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ABSTRACT 


This paper summarizes the Explicit Algebraic Stress Model in K-u form (EASM-ko) and in 
K-e form (EASM-ke) in the Reynolds-averaged Navier-Stokes code CFL3D. These models have 
been actively used over the last several years in CFL3D, and have undergone some minor 
modifications during that time. Details of the equations and method for coding the latest 
versions of the models are given, and numerous validation cases are presented. This paper 
serves as a validation archive for these models. 


IV 



1 INTRODUCTION 


A turbulence model must be employed in conjunction with the Reynolds-averaged Navier- 
Stokes (RANS) equations in order to close the equations. Over the last 30 years, many different 
turbulence models have been developed for use with RANS. Because there are so many to choose 
from, it is sometimes difficult for a CFD user to decide which model to use for a given flow 
situation. 

It is therefore very important that new (as well as existing) models be thoroughly validated 
in all the CFD codes into which they are implemented. This validation process should consist 
of applications to a suite of as many test cases as possible, which, taken together, thoroughly 
exercise the capabilities of the turbulence model. 

Most turbulence models for use with RANS today are linear eddy viscosity models, which 
assume a Boussinesq relationship between the turbulent stresses and mean strain rate tensor 
through the use of an isotropic eddy viscosity. Recently, however, nonlinear eddy viscosity 
models have been gaining widespread attention. This class of models assumes a higher-order 
tensor representation involving either powers of the mean velocity gradient or combinations of 
the mean strain rate and rotation rate tensors. 

One of the advantages of nonlinear eddy viscosity models over linear eddy viscosity models 
is that nonlinear models can predict differences in the turbulent normal stresses. Although this 
deficiency in linear models is not generally considered important for most external aerodynamic 
problems of interest, it has been shown to be crucial for capturing secondary motion in a corner 
flow, which is driven by the gradient of the turbulent normal stresses (see Abdol-Hamid et 
al.[l]). 

Explicit algebraic stress models (EASM) belong to the class of nonlinear eddy viscosity 
models. However, unlike some nonlinear models which determine their expansion coefficients 
through calibration with experimental or numerical data, EASMs obtain their expansion coef- 
ficients through their rigorous relationship with their “parent” full differential Reynolds stress 
equations. See Gatski and Rumsey [2] for details. 

The EASM models, when originally conceived [3], were solved by assuming a fixed (equilib- 
rium) value of V /e when determining the parameter g = (Ci/2 + V/e - l) -1 . This “fixed-#” 
EASM model was subsequently improved [4, 5, 6] by allowing V/e to vary. It is only this 
latter method, sometimes referred to as “variable-#” EASM, that is discussed in this report. 
The “fixed-#” EASM method has also been coded into CFL3D, but it is generally no longer 
recommended and will not be discussed. Other applications of EASM (using variable-#) can be 
found in Carlson [7] and Carlson et al.[8]. 

The purpose of this paper is primarily as a validation archive for the EASM-ko and EASM-ke 
models as coded in CFL3D. Thus, engineers who attempt to code these turbulence models for 
themselves in their own CFD codes can use this paper as a reference for the results they should 
expect to get with these models. This paper is not intended to advocate EASM, or to point 
out cases for which EASM yields improved (or worse) results than more conventional models. 
Therefore, no comparisons are presented between EASM and other models. That exercise is 
left to the reader. Here, the two versions of EASM-ko and EASM-ke are compared only with 
each other (and with experiment or theory). 

The next two sections describe the formulas for the EASM-ko and EASM-ke models, includ- 
ing details regarding the numerical method employed to solve them in CFL3D. Following that, 
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the models are applied to several test cases. Then, conclusions are drawn in the last section. 


2 OVERVIEW OF CFL3D 


The computer code CFL3D [9] solves the three-dimensional, time-dependent, Reynolds averaged 
compressible Navier-Stokes equations with an upwind finite-volume formulation (it can also be 
exercised in two-dimensional mode of operation for 2-D cases). It can solve flows over multiple- 
zone grids that are connected in a one-to-one, patched, or overset manner, and can employ 
grid sequencing, multigrid, and local time stepping when accelerating convergence to steady 
state. Upwind-biased spatial differencing is used for the inviscid terms, and flux limiting is used 
to obtain smooth solutions in the vicinity of shock waves, when present. Viscous terms are 
centrally differenced, and cross-diffusion terms are neglected. For very low Mach number flows, 
preconditioning [ 10 ] is used to insure convergence and accuracy of the solutions. 

The CFL3D code is advanced in time with an implicit approximate factorization method. 
The implicit derivatives are written as spatially first-order accurate, which results in block 
tridiagonal inversions for each sweep. However, for solutions that utilize Roe flux-difference 
splitting [ 11 ], the block tridiagonal inversions are further simplified using a diagonal algorithm 
with a spectral radius scaling of the viscous terms. All solutions in this paper use Roe. 

Following Wilcox [12], Reynolds averaging can be used with the Navier-Stokes equations in 
Favre variables [13] to account for turbulent fluctuations. The resulting equations of motion 
can be written using the summation convention as follows. The full Navier-Stokes equations are 
shown here, but in CFL3D they are solved using the thin-layer approximation in pre-selected 
coordinate direction (s). The equations are given here in terms of dimensional quantities. Al- 
though not shown, the equations are subsequently nondimensionalized and solved in generalized 
coordinates (see Krist et al.[9]). 


dp 

dt 


d . . 

= 0 


» , . 
«<'"“> + 


d dp 

(■ P u i u ' ) - 


dx 


dxi 


drjj 

dxj 


9 0 0 

—Sp e ) + - — (~ u i T a ~ v + tfj). 


dt 


dxi 


(1) 

( 2 ) 

(3) 


where 


p = (7 - 1 ) pE- ^p(u 2 + v 2 + w 2 ) - pK 

(4) 

E = e + ^{u 2 + v 2 + w 2 )+K 

(5) 

H = E + p/p 

( 6 ) 

• 1 / M , Pt\ da 2 

* 7 -lVPr Pr J dxj 

(7) 

a 2 =™. 

P 

(8) 


2 



where ok is defined in the next section. Define: 
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The stress term r,j is composed of a laminar and a turbulent component as 


T a = T ij + T ij> 

where 


( 10 ) 

( 11 ) 


(12) 



For linear eddy-viscosity models using the Boussinesq eddy viscosity assumption: 

rj = 3 pKS {j - 2 p* t (Sij - - |^,j) ( i4 ) 

where p* is the turbulent eddy viscosity computed by a linear turbulence model. For the 
nonlinear eddy-viscosity models EASM-ko and EASM-ke: 

T Ii = 3 pKSij - 2/4 (Sij - 3 

+[a 2 a 4 (SikWkj - W ik S k j) 

—2a3a 4 (SikSkj — -SkiSikSij)]) ( 15 ) 

where the terms a 2 , a 3, and a 4 are defined in the next section. 

In CFL3D, the kinetic energy of the fluctuating turbulent field K is ignored in the definition 
of p and E in equations (4) and (5). Furthermore, for all linear eddy- viscosity models, the term 
involving K is ignored in the definition of rfj in equation (14), and the tpj term in the energy 
equation (3) is assumed to be zero. These approximations and assumptions have relatively little 
effect for most flows at low or transonic Mach numbers, but could have significant impact for 
hypersonic flows. From a coding perspective, for linear eddy-viscosity models, the end result is 
that the laminar Navier-Stokes equations are identical to the turbulent Navier-Stokes equations 
with the exception that p is replaced by p + p* and p / Pr is replaced by p/Pr + p* /Pr t (in these 
equations Pr is taken as 0.72 and Prt is taken as 0.9). 

For the nonlinear models, the term involving K is included in the definition of tJ- in equation 
(15), and the rpj term in the energy equation (3) is modeled via equation (9). However, K is 
still ignored in the definition of p and E in equations (4) and (5). Therefore, like the linear 
turbulence models, strictly speaking the nonlinear turbulence models are not applicable for 



hypersonic flows either, as coded in CFL3D. From a coding perspective, the nonlinear models 
are similar to the linear models except that additional terms beside p* and p*/ Pr< must be 
added to the right-hand sides of the momentum and energy equations (2) and (3) when going 
from laminar to turbulent Navier-Stokes. 

All the turbulence models in CFL3D, including EASM-ko and EASM-ke, are solved uncou- 
pled from the mean flow equations using implicit approximate factorization. Their advective 
terms are solved using first-order upwind differencing. Details for obtaining the p* term for the 
EASM-ko and EASM-ke models are given in the next section. 


3 THE EASM EXPLICIT ALGEBRAIC STRESS MODELS 

For background and summary on the derivation of EASM, see Gatski and Speziale [3], Jongen 
and Gatski [14], Rumsey et al.[6], and Rumsey and Gatski [15]. The following section bypasses 
the derivation and jumps directly to the final resulting equations for the EASM-ko and EASM-ke 
models, and also outlines the method of solution in CFL3D. 

The kinematic eddy viscosity p* is given by 


p; = ti/p = C;KT=-Ka 1, 


(16) 


with r = l/u> for EASM-ko and r = K/e for EASM-ke. Thus, otifr is equivalent to -C*, 
which is variable (as opposed to constant as in most linear two-equation models). The nominal 
level for C* in a zero-pressure-gradient log layer is approximately 0.09. The value of ot\/r is 
obtained from the solution to the following cubic equation at each point in the flow field: 
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The correct root to choose from this equation is the root with the lowest real part [16]. 
Also, the degenerate case when r) 2 — >■ 0 must be avoided. See Rumsey and Gatski [15] for 
further details. In the current implementation, the resulting C* = — (cq/r) is limited by 
C* = max(-(ai/r), 0.0005). Other parameters are given by 
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(24) 

a2 = l -{ 2-c 4 ) 

(25) 

a 3 = - (2- C 3 ) 

(26) 

a 4 = 7j - 2y£ r) 2 r 2 t. 

(27) 

To = Cl/2 

(28) 


(29) 


and C e i = 1.44, C e2 = 1.83, C? = 3.4, C\ = 1.8, C 2 = 0.36, C 3 = 1.25, and C 4 = 0.4. 

The preceding implementation is exactly the same for EASM-ko or EASM-ke. The two 
methods differ in the two equations solved (K-u> vs. K-e). For EASM-ko, the explicit tensor 
representation for r tJ is coupled with the following K-u> two-equation model: 


DK 

Dt 


— V- fp'Ku + 


3 

dx k 



dlC 

dxk. 
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where 



( 31 ) 

(32) 


and ok = 1, <r w = k 2 /[^/UK/} — 7 )], k — 0.41, 7 = 0.53, /? = 0.83, and C ^ = 0.0895. Note that 
for 2-D incompressible flows, V = 2v*t] 2 is exact. In CFL3D, the user has the option of using 
the exact or the approximate V term. We have found there to be very little difference for a 
wide variety of subsonic and transonic 2-D aerodynamic-type flows, and use of the approximate 
term is more robust during transient stages of the computation in some cases. Except where 
otherwise noted, all the results in this paper were obtained using the approximate term. 


Also, it should be noted that the values of an and 7 in this model are different than reported 
in Rumsey and Gatski [15]. They were changed recently to improve the model’s capability for 
jet-type flows (see Georgiadis et al.[17]). The change was found to have relatively small impact 
for wall-bounded flows in general. In the current implementation, V in the A'-equation is limited 
to be less than 20 times the destruction term fp*Ku. The function , taken from from Wilcox 
[12], is given by 


fp, - 1 when Xk < 0 


l + 680 xg 

1 + 400x| 


when XJt > 0 


_ cl ok du 

* fc u ) 3 dxj dxj ’ 


(33) 
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where the C'l term in the formula for Xfc is necessary because ui in the current model does not 
“absorb” C ^ as in Wilcox’s model. 
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For the EASM-ke two-equation model: 


De 

Dt 


DK 

Dt 


= V -£ + 


dlC 
dx k . 


'.h%) 

= C„|p - f.cj- + £- k [(« + £) 


d 

dx k 

2 


de ‘ 
dx k \ ’ 


(36) 
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where V is given by equation (32) and is limited in the A'-equation to be less than 20 times 
the destruction term e. Also, f e = [ 1 — exp(-i?e/f/10.8)], Ren = A' 1 / 2 d/v , ctr = 1.0, cr e — 
k 2 — C e i)], C e i = 1*44, C e 2 = 1.83, C M — 0.0885, and d is the distance to the nearest 
wall. Additional wall damping functions (such as / M , to achieve expected asymptotic behavior of 
the turbulence quantities very near the wall) are not employed in the current EASM-ke model. 
Note that in equations (36) and (37), the diffusion terms are modeled using an equilibrium eddy 
viscosity Ut = C tl K 2 /e , where the constant C M = 0.0885 for this model. This is different than in 
the EASM-ko model, which uses the actual eddy viscosity v* (with variable C*) in its modeled 
diffusion terms. The diffusion terms for both EASM-ko and EASM-ke are approximate models 
in any case; see, for example, Warsi [18]. 


In CFL3D, the turbulence equations are nondimensionalized by the same reference variables 
as the Navier-Stokes equations: a re f, /> re f, Pref> and L re f, where a is speed of sound, p is density, p 
is molecular viscosity, and L is length scale corresponding to unit one of the computational grid. 
Thus, denoting nondimensional variables with a tilde: K — A/(a 2 ef ), Cj = (/i re fa;)/(p re fa 2 ef ), 

and E = (/^ref-) / (Pref^ref ) • 

The turbulent boundary conditions applied at solid walls are K w = 0, £ w = 2v w (d\/K /dn)l„ 
and u> w = 10(6i/u,)/[/?(An) 2 ], where An is the distance to the first cell center away from the 
wall. The boundary condition for u w is from Menter [19]. This boundary condition simulates 
the analytical behavior of u near solid walls without the need for specifying the solution at 
interior points. At farfield inflow boundaries, the following nondimensional values are assigned 
as freestream reference values in CFL3D: K = 9 x 10 -9 and u = 9 X 10 -8 for EASM-ko; 
and K = 1 x 10 -9 and § = 1 x 10 -17 for EASM-ke. The end result for both models is that 
nondimensional fit = Pt/Pref « 0.009 at farfield inflow boundaries (based on C* = 0.09, which is 
only approximate because C* is variable in these models). At outflow boundaries, the turbulent 
quantities are extrapolated from the interior of the computational domain. 

Some additional details concerning the method for solving the two-equation turbulence 
models are now given. First, the following is an example showing the transformation from 
Cartesian coordinates to generalized coordinates: 


where 


U .<^ =U M +V dk dk 

3 dxj d$ + dr) + dC 


U = £ x u + £ y v + £ z w + Zt 


(38) 


(39) 


V = Tj x U + T)yV + T) Z W + rjt 


(40) 


W = ( X U + C yV + C Z W + Ct 


(41) 
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and £, 77 , and £ are the three generalized coordinate directions. When transforming into gener- 
alized coordinates, cross-derivative terms are neglected in the diffusion-type terms that appear 
in the turbulence equations. 

The two turbulence equations are solved decoupled using implicit approximate factorization. 
Each sweep requires the solution of a scalar tridiagonal matrix. The production terms are 
treated explicitly, lagged in time, while the destruction and diffusion terms are treated implicitly. 
It should be noted here that in CFL 3 D, the K, w, and e terms are not allowed to go negative 
during the solution procedure. If the update yields a negative result, it is instead limited to 
a very small positive number, and a counter keeps track of the number of times this occurs 
in the flowfield as an indicator of non-convergence in the solution. Furthermore, rji, T 22 , and 
733 are limited to be positive (this is the realizability constraint). In EASM-ke, t) 2 t 2 and 
TZ 2 t] 2 t 2 are each limited to be less than 1200 in magnitude. In both EASM-ko and EASM-ke, 
where uor£ appears in the denominator of equation (15) (i.e., in the nonlinear terms added 
to the Navier-Stokes equations), these values are limited to be greater than or equal to their 
freestream reference levels in order to avoid dividing by unrealistically small numbers during 
the convergence process. 

4 VALIDATION TEST CASES 

This section contains a detailed account of a series of validation test cases to which the EASM- 
ko and EASM-ke models have been applied. The first 7 cases are also “standard” CFL3D test 
cases, described and available on the CFL 3 D website [20]. (In some instances, the case on the 
website has been set up for a different turbulence model other than EASM-ko or EASM-ke, 
and some adjustments to CFL number are necessary in the provided input file in order to get 
the cases to run.) The last 2 cases are not currently available on the website. For the first case 
(flat plate), grid sensitivity studies are described in this report. All of the subsequent cases use 
only a single grid size, with no parametric variations performed. However, many of the grids 
used in these validations were found to be fine enough to capture the flow features of interest 
in previous grid sensitivity studies performed by the authors or by other researchers. In each 
case, the boundary conditions are described in detail; unless otherwise noted, at all solid walls 
a no-slip, adiabatic wall boundary condition is used. At these walls, the turbulence boundary 
conditions are those described in section 3. For the mean flow variables, velocity components 
at the wall are set to zero, dT w /dn = 0, and p w is determined using linear extrapolation from 
the interior of the domain. Note that the EASM-ke model sometimes has difficulty establishing 
turbulence when started from freestream or weak turbulence initial conditions. Therefore, for 
all the results in this paper, EASM-ke was initialized from the corresponding EASM-ko solution. 

4.1 Flat Plate 

The flat plate case (zero pressure gradient) was performed at M = 0.3, Re = 6 X 10 6 per unit 
length of the grid. The default grid size for this case was 65 X 97. Fig. 1 shows a picture of 
the grid. The viscous, adiabatic wall was unit one in length. Symmetry boundary conditions 
were imposed on a region of length x = 0.333 in front of the plate. Grid domain height 
was approximately 1 . 0 , minimum normal grid spacing at the wall was 1 x 10 -6 , and the grid 
was stretched at a rate of 1.18 until the vertical spacing first exceeded the horizontal spacing. 
Horizontal grid spacing was constant at Ax = 0.0208333. 
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At the inflow boundary, the total pressure and total temperature were set via pr/Pref = 
1.02828 and Tj/T re i = 1.0080, and a 1-D characteristic boundary condition was used to deter- 
mine the other flow variables [9]. At the outflow boundary, linear extrapolation was employed. 
At the top boundary, a farfield Riemann-invariant 1-D characteristic boundary condition was 
employed. 



Figure 1: Flat plate grid, 65 x 97. 

Surface skin friction coefficient results using EASM-ko are shown in Fig. 2 for a grid density 
study. A total of four grids were used: each had 65 points in the streamwise direction, and the 
number of points in the normal direction was successively halved from 193 to 97 to 49 to 25. 
Thus, the default grid size (65 X 97) was included as well as one level finer and two levels coarser. 
Also shown in the figure are symbols giving three different theoretical curves from White [21], 
indicative of the type of spread that might be expected from experimental data of a turbulent 
flat plate boundary layer. As the grid is refined, the results approach the range defined by the 
theoretical curves. Looking at one particular x-location on the plate (x = 0.8125), and using 
second-order Richardson extrapolation on the finest two grids, one finds that the 65 x 193 grid 
yields a skin friction value 0.5% in error from the result on an infinitely-refined grid at this 
location. The default 65 x 97 grid is 2.0% in error, the 65 x 49 grid is 4.7% in error, and the 
65 x 25 grid is 10.6% in error. 

The same grid study was also conducted for EASM-ke. Results are shown in Fig. 3. Unlike 
EASM-ko, with EASM-ke the skin friction decreases with increasing grid density. In this case 
the 65 X 193 grid yields a skin friction value 0.1% in error from the result on an infinitely-refined 
grid at x = 0.8125. The default 65 x 97 grid is 0.6% in error, the 65 x 49 grid is 2.9% in error, 
and the 65 X 25 grid is 9.0% in error. Notice that the infinitely-refined skin friction levels 
approached by EASM-ke are slightly lower than those using EASM-ko. 
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Figure 2: Effect of grid size normal to the wall on surface skin friction, EASM-ko. 
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Figure 3: Effect of grid size normal to the wall on surface skin friction, EASM-ke. 
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Figs. 4 and 5 show the effect of varying the minimum normal spacing at the wall using 
EASM-ko and EASM-ke, while keeping grid size constant (at 65x97). The “standard” minimum 
spacing of 1 X 10 -6 corresponds with an average y+ value for the first cell-center off the wall 
of approximately 0.2. The minimum normal spacing was decreased to 1 x 10 -7 (average y+ of 
0.02) as well as increased to 2.2 x 10 -6 (average y+ of 0.5), 5 X 10" 6 (average y-|- of 1.1), and 
1 x 10 -5 (average y+ of 2.2). In each case the stretching rate was varied to keep the grid height 
approximately 1.0. As the average y+ level decreases, the skin friction increases for EASM-ko 
and decreases for EASM-ke. Using the result at x = 0.8125, the difference between the smallest 
y+ and the largest y+ tested is approximately 5% for EASM-ko and 2.6% for EASM-ke. The 
difference between using y+ of 0.02 and 0.2 is only 0.3% for EASM-ko and less than 0.1% for 
EASM-ke. 

A plot of u+ vs. log(y+) is shown for the two models in Fig. 6 using the 65 x 97 grid, along 
with Spalding theory from White [21]. Results agree well with the law-of-the-wall. 



X 


Figure 4: Effect of average minimum y+ on surface skin friction, EASM-ko, 65 x 97 grid. 


4.2 Back Step 

The back step case compares against experimental data from Driver and Seegmiller [22], In the 
experiment, several different upper wall angles were tested, but in the present computation, the 
upper wall was straight (0 deg). The test conditions were M = 0.128, Re = 37573 per unit step 
height H. In the grid for this configuration, the lower wall starts out at a height of 1.0, then 
drops to a height of 0.0 at x/H = 0. The upper wall is at a height of 9.0. The grid extends 
from x/H = -4.0 upstream to x/H = 35 downstream. 
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X 


Figure 5: Effect of average minimum y+ on surface skin friction, EASM-ke, 65 x 97 grid 
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Figure 6: Velocity profiles for both models in wall variables at x = 0.8, 65 x 97 grid, y+=0.2. 
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The grid itself is made up of two zones: 25 x 65 and 129 x 113. These two zones abut along 
x = 0 in a patched-grid (non 1-to-l) fashion. Interpolation is used in CFL3D at the patched 
interface to provide communication between the zones. The minimum normal spacing at the 
walls varies within the range of approximately 0.006-0.007, which yields an average y+ value at 
the first cell center off the walls of between 0.6 and 0.9. The only exception to this is the back 
face of the step itself, which has a normal spacing of 0.1. The streamwise grid spacing varies. 
It is clustered near the step, and stretched as it approaches the upstream and downstream 
boundaries. Two pictures of the grid are given in Figs. 7 and 8. 

The lower wall, upper wall, and back of the step used viscous, adiabatic boundary conditions. 
At the upstream boundary, the density, velocity, and turbulence quantities were set to closely 
match experiment (see Rumsey et al.[6] for details concerning this procedure). The pressure 
at the inflow was extrapolated from the interior of the domain. At the downstream boundary, 
the pressure was set to p/p re { = 1.00149 and all other quantities were extrapolated from the 
interior of the domain. 



Figure 7: Backstep grid, 25 x 65 and 129 x 113. 



Figure 8: Close-up of the backstep grid near the interface between zones. 


12 



Lower wall pressure coefficient and skin friction coefficient for both EASM-ko and EASM-ke 
are given in Figs. 9 and 10. A series of velocity profiles is given in Figs. 11, 12, and 13, and 
turbulent shear stress is given in Figs. 14, 15, and 16. 

Reattachment length in the experiment was approximately in the range of x = 6 — 6.3. 
Using EASM-ko, the computed reattachment length was x = 6.30, whereas for EASM-ke it was 
x ~ 5.86. Both of these are in very good agreement with the experiment, considering that many 
turbulence models have been shown to seriously underpredict reattachment length by 30% or 
more [22], The EASM-ko skin friction results agree extremely well with experiment in general, 
but EASM-ke underpredicts the levels in the separated region and overpredicts the levels after 
reattachment. This poorer prediction of EASM-ke is consistent with skin friction results that 
others have found for this case using different k-e models (see, e.g., Menter [19]). Velocity 
profiles for both models are similar. Both do a reasonably good job in the separated region, 
but then do not predict as rapid a recovery after reattachment as experiment. This too-slow 
recovery from separation is a chronic problem with almost all RANS turbulence models in use 
today [23]. Turbulent shear stress is predicted well by both models up through x = 2.5, but 
then the peak shear is underpredicted (and is also farther from the wall) through approximately 
x = 5.5. Downstream of this, the levels and shapes again agree fairly well with experiment. 



Figure 9: Backstep surface pressure coefficients. 


4.3 Transonic Diffuser 

This case models the strong-shock diffuser experiment of Sajben and Kroutil [24]. Fig. 17 shows 
the 81 X 51 grid, which was obtained from the NPARC Alliance Validation Archive website [25]. 
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Figure 16: Backstep turbulent shear stress profiles, set 3. 

This configuration had an entrance to throat area ratio of 1.4 and an exit to throat area ratio 
of 1.5. This flow is characterized by the ratio, R, of exit static to inflow total pressure. For 
the strong-shock case R is 0.72. In the strong shock case, there is separated flow on the upper 
wall downstream of the shock wave [25]. The grid had a minimum spacing at the walls of 
approximately 0.0001, yielding an average minimum y+ value at the first cell center off the wall 
of approximately 2.5. 

When computing this case in CFL3D, the lower and upper walls were modeled as viscous, 
adiabatic walls. At the inflow plane, the total pressure divided by reference pressure was 1.894 
and total temperature divided by reference temperature was 1.201. At the outflow plane, the 
static pressure divided by reference pressure was 1.364. Reference conditions (at the throat) 
were: M=1.0 and Re=937,000 per throat width. (In the grid, the throat width was unit 1.) 

Computed wall pressures are shown compared to experiment in Figs. 18 and 19 for the lower 
and upper walls, respectively. There is a slight difference in predicted shock position, with 
EASM-ko downstream of EASM-ke, but both models agree well with experiment in general. 
Velocity profiles at four x-locations are shown in Fig. 20. Both EASM-ko and EASM-ke yield 
similar results in good agreement with experiment. Both models also correctly predict separated 
flow on the upper wall at the two stations x/H=2.662 and 4.611, although in both cases the 
maximum reverse flow velocity magnitude is lower than in the experiment. 
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Figure 19: Transonic diffuser upper wall pressures. 

4.4 NACA 4412 Airfoil 

This case models the separated flow over a NACA 4412 airfoil, and compares with the exper- 
iment of Coles and Wadcock [26]. This flow was modeled computationally using a freestream 
Mach number of M = 0.2, Reynolds number Re = 1.52 X 10 6 based on chord, and angle of 
attack of a = 13.87°. This flow separates on the rear part of the airfoil upper surface at these 
conditions. Fig. 21 shows a picture of the 257 X 81 grid used. The grid had a minimum spacing 
at the wall of 2 X 10 -5 chords, a farfield extent of approximately 17 - 20 chords, and a total of 
177 points on the airfoil surface. The average minimum y+ value of the first cell center off the 
wall was less than 2 for this grid. 

Surface pressure coefficients are shown in Fig. 22 and upper surface velocity profiles are 
shown in Fig. 23. Both turbulence models predict similar surface pressures, which are in good 
agreement with experiment everywhere except near the upper surface trailing edge. For velocity 
profiles, EASM-ko yields results in excellent agreement with experiment, whereas EASM-ke 
predicts velocity levels that are too low (the profile has greater inflection). However, both 
models predict separation at the same location of x/c = 0.830. 

Part of the reason for the poorer prediction of EASM-ke may be due to the fact that the K-e 
formulation in general does not handle wall-bounded adverse pressure gradients well, because 
of shortcomings in the e equation. See, for example, Rumsey and Gatski [15], Wilcox [12], 
and Rodi and Scheuerer [27]. Rodi and Scheuerer concluded that the generation term of the 
s equation has to be increased above its usual level. Fig. 24 shows a plot of the predicted 
velocities in wall variables as compared to Spalding theory at x/c = 0.5 on the upper surface 
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of the NACA 4412 airfoil. EASM-ke did not obtain the correct slope of the log layer in this 
adverse pressure gradient part of the flow, whereas EASM-ko did. This same problem was also 
shown and discussed in Rumsey and Gatski [15] and Nagano and Tagawa [28]. 



Figure 21: NACA 4412 grid, 257 X 81. 


4.5 Ejector Nozzle 

The ejector nozzle case models a subsonic 2-D jet flow that entrains and mixes with a secondary 
outer flow. The experimental data are from Gilbert and Hill [29]. See also Georgiadis et al.[30]. 

This case has nozzle total pressure divided by atmospheric static pressure of 2.44, with 
total temperature = 644 R. The secondary flow has a total pressure equal to atmospheric 
and a temperature of 550 R. At the outflow used by this grid, the static pressure divided 
by the atmospheric static pressure was taken to be 0.9131. However, in CFL3D, a reference 
condition was chosen that corresponded with a non-zero Mach number: in this case it was 
taken to be M=0.22, which is the approximate Mach number at the secondary flow inlet in 
this grid. (Corresponding Reynolds number was taken to be approximately 1.64 million per 
foot.) Using M=0.22 and isentropic relations, total pressure divided by reference pressure is 
1.0343 and total temperature divided by reference temperature is 1.0097. These numbers are 
used as inflow conditions for the secondary flow. Also as a result of choosing this reference 
condition, static pressure divided by reference pressure at the outflow is 0.9444. At the nozzle 
inlet, total pressure divided by reference pressure is 2.5237 and total temperature divided by 
reference temperature is 1.1822. 

The grid consists of three zones with the following dimensions: 31 x 41 for primary nozzle, 
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Figure 22: Surface pressure coefficient for NACA 4412. 

31 X 71 for secondary inlet, and 101 X 121 for downstream. The grid is shown in Fig. 25. 
Minimum grid normal spacing at the walls yielded an average y+ level of approximately 0.7. 

Fig. 26 shows velocity profiles at four stations downstream using EASM-ko and EASM- 
ke. Both models overpredict the centerline velocity at the upstream stations and somewhat 
underpredict the jet width in general, but qualitatively the overall agreement is good (it is 
similar to results on the same grid using different turbulence models and a different code in 
Georgiadis et al.[30]). 

4.6 Axisymmetric Bump 

The axisymmetric bump case has served as a standard test case for transonic separated flow for 
many years. The experimental data are from Bachalo and Johnson [31], and have been used, 
for example, to validate Menter’s widely-used SST turbulence model [32]. 

In this case, a turbulent boundary layer develops in the axial direction over an axisymmetric 
circular bump (smoothed at its leading and trailing edges). Wall effects in the experiment are 
minimal because the body does not generate any lift. For the particular case studied here, the 
Mach number is 0.875 and Reynolds number based on the length of the bump is 2.66 X 10 6 . At 
these conditions, there is a shock wave and the flow separates on the rear part of the bump, 
then reattaches a short distance downstream. 

This case is run in “3-D mode” in CFL3D on a grid with two planes separated by 1.0 degrees, 
and with periodic boundary conditions used on each of the two circumferential planes. Fig. 27 
shows a picture of one plane of the 181 x 101 x 2 grid. The grid had an approximate minimum 
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Figure 24: Velocity profiles in wall variables for upper surface of NACA 4412 at x/c — 0.5. 


secondary inlet 



Figure 25: Ejector nozzle grid. 
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spacing at the wall of 1.5 x 10~ 5 chords. Its farfield extent was 3.2 chords in front of the 
bump leading edge, 3.4 chords downstream of the bump trailing edge, and 4.1 chords domain 
height. There were 85 axial points on the bump part of the wall surface. The average minimum 
y+ value of the first cell center off the wall was approximately 1.3 for this grid. Farfield-type 
boundary conditions were used on the upstream, downstream, and top faces of the grid. 

Fig. 28 shows surface pressure coefficients on the body. Both EASM-ko and EASM-ke 
capture the shock position very accurately, as well as the shape of the pressure variation down- 
stream. EASM-ke predicts the shock location slightly further forward. Velocity profiles and 
turbulent shear stress profiles at seven stations are shown in Figs. 29 and 30, respectively. At 
the station xjc — 0.688 upstream of separation, velocity profiles using EASM-ke agree with 
the experimental data somewhat better than EASM-ko, but EASM-ko recovers better than 
EASM-ke at the last three stations, downstream of reattachment. Both models yield excellent 
velocity profile results in the separated region. Also, both models predict turbulent shear stress 
profiles in reasonable agreement with experiment, with EASM-ke tending to give slightly lower 
peak levels in magnitude in the separated region. 

In section 3, the use of an approximate production term in the EASM-ko and EASM-ke 
models was discussed, and it was asserted that use of the approximate term makes very little 
difference for a wide variety of subsonic and transonic 2-D aerodynamic- type flows. In fig. 31, 
shear stress profiles at two stations are shown using the two models with both the approximate 
V term as well as its exact form. Results are nearly identical. Pressure coefficients and velocity 
profiles (not shown) are also nearly identical. 



Figure 27: Axisymmetric bump grid, 181 x 101 x 2. 


26 



1 

x/c 


Figure 28: Surface pressure coefficients for axisymmetric bump. 

4.7 ONERA M6 Wing 

The ONERA M6 wing is a 3-D case often used to validate codes and turbulence models. The 
experimental data are from Schmitt and Charpin [33]. 

The particular case studied here has a Mach number of 0.84 and Reynolds number based on 
mean aerodynamic chord of 11.7 x 10 6 . The angle of attack is 3.06°. At these conditions, there 
are two shock waves across the inboard part of the wing upper surface, which merge to form a 
single shock wave at the outboard part of the wing. There is no significant flow separation for 
this case. 

The grid used for this case is a C-0 type grid of size 289 x 65 x 49. Part of the grid (with 
every other grid point removed) is shown in Fig. 32. The grid is scaled so that the chord of the 
wing varies from 0.6737 at the root to 0.3789 at the tip. In these units, the Reynolds number 
per unit length is 21.66 X 10 6 . The grid extends to approximately -6.4 unit lengths upstream, 
7.4 unit lengths downstream, 6.4 unit lengths high, and 7.4 unit lengths to the side. Symmetry 
conditions are applied along the center plane, and farfield conditions are applied at the farfield 
boundaries. There are 225 streamwise points and 49 spanwise points on the wing surface itself, 
for a total of 11025 surface points. 

The average minimum normal spacing at the wall is approximately 6 x 10 -6 unit lengths, 
which yields an average minimum y+ value of the first cell center off the wall of approximately 
3.2. This minimum spacing is larger than one would usually want to use for a turbulent flow 
computation. Based on the flat plate results in section 4.1, such a large wall spacing can yield 
skin friction levels that are more than 5% in error. However, it was still considered important 
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Figure 31: Effect of approximate production term on turbulent shear stress profiles for axisym- 
metric bump. 

to include this case as a part of this validation archive because it demonstrates the capabilities 
of these models on a less-than-optimal grid (which is not an uncommon occurrence in many 
situations involving 3-D configurations). 

Fig. 33 shows surface pressure coefficients at six span stations for the two models. Both do 
an excellent job capturing the pressure levels, including shock location. 

4.8 2-D Airfoil Wake 

This test case investigates the near-field development and decay of the 2-D wake of an airfoil 
at low Mach number. The experimental data are from Nakayama [34]. 

The case studied here uses the “Model A” airfoil, which is a 10% thick conventional airfoil. 
It is computed at a Mach number of 0.1 and Reynolds number based on chord of 1.2 x 10 6 
(the chordlength is 24 inches). The angle of attack is 0°. At these conditions, there is no flow 
separation on the airfoil. 

The grid used for this case is a C-grid, but it is divided into two zones. The first zone is 
257 x 97 and consists of the part of the C-grid containing points on the airfoil surface (thus, 
there are 257 points on the airfoil surface). The second zone is 153 x 193 and consists of the part 
of the C-grid in the wake, with the top and bottom halves put together. The grid is shown in 
Fig. 34. The grid extent is approximately 20 chords, and the average minimum normal spacing 
at the wall is approximately 1 x 10 -5 chords, which yields an average minimum y-f value of 
the first cell center off the wall of approximately 0.5. In the wake, the minimum grid spacing 
occurs along the wake cut; it starts out with the same minimum spacing as at the airfoil, then 
spreads further downstream. For example, at x/c = 1.2, the minimum normal spacing on the 
centerline is approximately 1 x 10 -4 chords. At that station there are approximately 97 grid 
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Figure 32: Schematic of ONERA M6 grid, 289 x 65 x 49 (every other grid point shown). 

points between y = — 1 and 1 inches ( y/c = —0.0417 and 0.0417). 

Fig. 35 shows wake velocity profiles at seven streamwise stations for the two models. The 
EASM-ke has a slightly deeper maximum wake deficit than EASM-ko at stations downstream of 
x/c = 1.2, in better agreement with experiment, but EASM-ko yields a slightly better spreading 
rate in this case. Both models predict the wake position to be somewhat too low compared to 
experiment at the last two stations. 

Fig. 36 shows wake turbulent shear stress profiles. Overall, both models agree well with 
experiment. The EASM-ke predicts higher peak levels than EASM-ko at the first two stations, 
then predicts slightly lower levels downstream. Both models underpredict the peak levels seen 
in the experiment at most of the stations downstream of xjc — 1.2. 

4.9 Curved Duct in Zero Pressure Gradient 

The curved duct test case was the subject of an extensive study in Rumsey et al.[35] related 
to curvature effects, and included the application of curvature corrections both to EASM-ke as 
well as to a linear one-equation model. However, the subject of the present paper is solely a 
validation for the “standard” forms of EASM-ko and EASM-ke, so curvature-corrected results 
are not presented here. (Nonetheless, it was found in Rumsey et al. that the “standard” form of 
EASM performs reasonably well for this case because EASM still retains some of the invariance 
properties of the full differential stress model. The more exact curvature-corrected version 
yielded only modest improvements.) The experimental data for this case are from So and 
Mellor [36]. 
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Figure 34: Nakayama model A grid, 257 x 97 and 153 x 193. 

This case is an internal flow with inner-wall convex curvature. The outer wall is shaped 
to yield approximately a zero streamwise pressure gradient along the inner wall in the region 
of interest. The Mach number at the inlet is 0.063, and the Reynolds number is taken to be 
3.6417 x 10 4 per inch. 

The grid used for this case is 129 x 81, shown in Fig. 37. In this figure (and in those to 
follow), the coordinate s refers to distance measured along the inner wall, with s = 24 in. as 
the location of the inlet to the computational grid. A grid study was performed in Rumsey et 
al.[35], where it was concluded that this grid size is adequate. Slip- wall boundary conditions 
were applied at the outer wall in the CFD simulation. This boundary condition allowed the 
simulation to be run without the complication of having to contend with tangential air or bleed 
boundary conditions. The inner wall boundary condition was no-slip, adiabatic wall. The 
minimum normal spacing at the lower wall was 0.00015 in., which yields an average minimum 
y+ value of the first cell center off the wall of less than 0.5. The outflow boundary condition 
set pressure at p/p re f = 1-0 (where “ref’ refers to inlet conditions), and extrapolated all other 
quantities. At the upstream boundary, the density, velocity, and turbulence quantities were set 
to closely match experiment (see Rumsey et al.[6] for details concerning this procedure). The 
pressure at the inflow was extrapolated from the interior of the domain. 

Fig. 38 shows surface skin friction coefficient along the inner wall. Both EASM-ko and 
EASM-ke give reasonable agreement with experiment. Velocity profiles (referenced to inlet 
conditions and local boundary layer thickness) are shown in Fig. 39. These are predicted in 
good agreement with experiment. Turbulence quantities are given in Figs. 40, 41, and 42. Note 
that all turbulent shear and normal stresses are in the local body/normal coordinate system. 
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Figure 35: Wake velocity profiles for Nakayama model A. 
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Figure 36: Wake turbulent shear stress profiles for Nakayama model A. 
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Stresses in this frame are related to those in the Cartesian frame by the following relations: 

u'v' = ^(v'v'c — wV c )sin(20) + wV c cos(20) (42) 

£ 

u'u 1 = u'u' c cos 2 0 + uV c sin 2 0 + uV c sin(20) (43) 

v'v' = vV c cos 2 0 4- «V c sin 2 0 - uV c sin(20), (44) 

where the subscript c indicates Cartesian frame, and 0 is the angle that the body tangent 
vector makes with the z-axis. 

The u’u’ turbulent normal stress levels are predicted somewhat too large in magnitude by 
the EASM models in the curved region of the duct, but the turbulent shear stresses and the 
v’v’ normal stresses are predicted in good agreement with experiment. 



x, inches 

Figure 37: Curved duct grid, 129 x 81. 


5 CONCLUSIONS 

In conclusion, this paper summarized the Explicit Algebraic Stress Model in K-u form (EASM- 
ko) and in K-e form (EASM-ke) in the Reynolds-averaged Navier-Stokes code CFL3D. Details 
of the equations and method for coding the latest versions of the models were given, and 
numerous validation cases were presented. Except for the well-known problem of the standard 
form of the e equation being ill-suited for wall-bounded adverse pressure gradient flows, both 
of these models were shown to yield good results for a wide variety of different cases. These 
cases included flow fields with shock waves, curvature, and significant regions of separation. 
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Figure 41: Turbulent u’u’ normal stress profiles, referenced to inlet conditions. 
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Figure 42: Turbulent v’v’ normal stress profiles, referenced to inlet conditions. 
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